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Abstract 

Interstitials and vacancies in the Abrikosov phase of clean Type II superconductors 
are line imperfections, which cannot extend across macroscopic equilibrated samples at 
low temperatures. We argue that the entropy associated with line wandering nevertheless 
can cause these defects to proliferate at a sharp transition which will exist if this occurs 
below the temperature at which the crystal actually melts. Vortices are both entangled and 
crystalline in the resulting "supersolid" phase, which in a dual "boson" analog system is 
closely related to a two-dimensional quantum crystal of He'^ with interstitials or vacancies in 
its ground state. The supersolid must occur for B ^ Bx , where Bx is the decoupling field 
above which vortices begin to behave two-dimensionally. Numerical calculations show that 
interstitials, rather than vacancies, are the preferred defect for B ^ (po/X\, and allow us to 
estimate whether proliferation also occurs for B ^ Bx- The implications of the supersolid 
phase for transport measurements, dislocation configurations and neutron diffraction are 
discussed. 
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1. INTRODUCTION 



Fluctuations, especially in high-temperature superconductors, play a prominent role 
in determining vortex configurations in Type II materials in an external field [1]. It now 
appears, for example, that clean single crystal samples of YBCO (in the absence of twin 
boundary pinning) melt via a first-order phase transition [2] at a temperature well 
below the upper critical field line B.c2{T) predicted by mean field theory [3-5]. Point 
disorder, in the form of oxygen vacancies, does not seem to affect this phase transition 
strongly in YBCO. It is quite possible that the disorder-induced translational correlation 
length [6] Ra greatly exceeds the vortex spacing for T < in the field range for which 
the transition is first order. Although in the presence of randomness the low temperature 
phase is not, strictly speaking, a solid, the thermodynamic properties at and near the 
phase transition should be similar to those in the absence of randomness. 

Consider the thermal excitations about the crystalline state on scales less than Ra- 
We assume for simplicity vortices which are perpendicular on average to the Cu02 planes, 
i.e. parallel to the ^-direction. The finite reduction of the Debye- Waller factor associated 
with the translational order parameter 

= exp [-i GiGj{ui{f±,z)uj{f±,z))] (1.1) 

by phonons is discussed in Ref. 7. Here, G is a reciprocal lattice vector and u{f±,z) is 
the displacement field of a flux lattice with vortices parallel on average to the z direction. 
Dislocation loops are a topologically distinct excitation which, when they proliferate at a 
melting transition, drive Pq{T) to zero and can lead to a hexatic fliix liquid with residual 
bond orientation order [8]. Isolated dislocation loops are far more constrained than their 
counterparts in crystals of point particles: dislocation loops in fact must lie in a plane 
spanned by their Burgers vector and the average field direction [9] — see Fig. la. 

Vacancies and interstitials differ even more dramatically from the analogous defects in 
crystals of point particles. The number of flux lines is conserved, which means that these 
defects are lines instead of points. The point- like nature of vacancies and interstitials in 
conventional crystals ensures that they are present in equilibrium at all finite temperatures 
for entropic reasons [10]. However, because such imperfections have an energy proportional 
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to their length in fliix crystals, they cannot extend completely across an equilibrated macro- 
scopic sample at low temperatures. A typical fluctuation at low temperatures might consist 
of the vacancy-interstitial pair shown in Fig. lb. Unlike the dislocation loop in Fig. la, 
this loop is not constrained to lie in a single plane. This configuration also provides a 
mechanism by which vortices near the loop may jog one lattice constant to the right as z 
varies. 

Defect loops of this type can have important dynamical consequences. In two dimen- 
sions, point-like vacancies and interstitials probably dominate the resistive properties of a 
weakly pinned lattice in a superconducting film [13, 14]. The same may well be true of 
vacancy-interstitial loops and lines in three dimensions as discussed in Section 5. 

Although the energy of vacancy and interstitial lines is proportional to their length, 
it is nevertheless possible for these defects to "proliferate" (i.e., to become infinitely long) 
at high temperatures for entropic reasons. Consider, for example, a vacancy wandering 
across a macroscopic sample of thickness L, as in Fig. 2. To estimate the free energy of 
this defect, we describe its trajectory along the z axis by a function fd{z) and write its 
partition function as a functional integral, 
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(1.2) 



Here Ed is the energy of an isolated defect and Sd is its tilt energy, defined in analogy 
with similar quantities for isolated fiux lines near Hd [7] and Ui{fd) is a periodic lattice 
potential with minima at the sites of the triangular Abrikosov crystal. Implicit in the path 
integral (1.2) is a length scale £z which is the average distance along z between hops of the 
vacancy from one lattice position to another. As a crude estimate of the path integral we 
replace it by exp(— e^L/T) 6^/^^, since the vacancy has six directions in which to hop on 
a triangular lattice. The free energy Fd = —T In Zd is thus 

Fd ^ SdL - ^ Lln6 , (1.3a) 

which becomes negative for T > Td, where 

rd = £d4/ln6. (1.3b) 

Above this temperature (provided the crystal does not melt first), vacancies (or intersti- 
tials) will proliferate in a crystalline phase. If these defects do not become strongly pinned 



3 



by point disorder, one might then expect a linear contribution to the resistivity due to 
defect motion within pinned bulk crystaUites, and a vortex glass transition [12] for weak 
pinning might then occur near T^. More precise estimates of the proliferation temperature 
will be represented later in this paper. A related phenomenon was suggested by Feigel'man 
et al. [13], who predicted the unbinding of "quarters" (i.e., quartets) of dislocations above 
the decoupling field Bx in highly anisotropic superconductors (see below). 

The phase in which defects such as interstitials and vacancies proliferate is in fact both 
crystalline and entangled. Regarded in light of the analogy between thermally excited flux 
lines and two-dimensional bosons [15], it represents a "supersolid" quantum crystal, in 
which vacancies and interstitials are incorporated into the ground state [15-20]. The 
possibility of a supersolid phase for flux lines in Type II superconductors was first noted 
on the basis of the boson analogy by M.P.A. Fisher and Lee [21]. Somewhat paradoxically, 
the "supersolid" is actually less superconducting for the BCS condensate electrons than a 
conventional vortex lattice phase. An alternative name is an "incommensurate solid" or a 
"vortex density wave." 

Figure 3 illustrates the connection between defect proliferation and the boson order 
parameter for flux lines [22] . This figure represents a low-temperature contribution to the 
order parameter correlation function 



where i(^{r±, z) and il)*{f'^, z') are destruction and creation operators for flux lines moving 
along the z axis, which plays the role of "time" for the "bosons." The composite operator 
in Eq. (1.4) creates an extra line at (r^,z') and destroys an existing line at {r^,z). The 
lowest energy configuration is then a line of vacancies (for z' > z, as in Fig. 3) or interstitials 
(for z' < z) connecting (r^, z') to {f±, z). At low temperatures, the defect line joining the 
head to the tail costs an energy of order /d(.s)s, where s = -^J {r± — r[_)'^ + {z — z')^ is 
the length of the string in direction s and fd{s) is the angle-dependent defect-free energy 
per unit length. The correlation function (1.4) then decays exponentially to zero, like 
exp(— /(/s/T), as |rx — — > oo with z and z' fixed. Above T^, fd changes sign, defects 
proliferate, and there is long-range order in G{f±, r^; z, z'), 



G{r^,f:^;z,z') = (V;(r-1, ^)V;*(fi, z')) , 



(1.4) 




lim G{f±,f^;z,z') = {i;{f±, z)){ij*{f^, z')) 



(1.5) 
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The physical meaning of Eq. (1.5) is that the free energy of the extra line segment 
from {f,z) to {r'jz') in the limit of large separations remains finite. The quantity — InG 
measures the free energy of a monopole anti-monopole pair at a (r, z) and {r', z') respec- 
tively. In the conventional solid phase the magnetic monopoles are confined by a linear 
potential, while when the defects proliferate, the monopoles are unconfined as in a normal 
metal with finite total free energy cost. Entanglement of vortex lines in the crystal will 
be catalyzed by the proliferation of vacancies and interstitials, since these allow fluxons 
to easily move perpendicular to the z axis. Of course, entanglement will also arise if the 
crystal melts directly into a flux liquid. In either case the boson order parameter becomes 
nonzero, 

MT) = mr±, z)) ^ . (1.6) 

Note that Eq. (1.6) gives a precise meaning to the concept of "entanglement" as 
used here and elsewhere to describe "superfluid" phases within the boson analogy. Strictly 
speaking, the concept of "tangled vortex lines" does not distinguish sharply between low 
and high temperature phases for lines in an infinitely thick sample. Indeed, in all phases 
discussed here — solid, "supersolid" and fluid — a vortex line will wander as a function of 
2; as a random walker on sufficiently large scales. Although this wandering will be more 
pronounced in the higher temperature phases, it does not uniquely distinguish them from 
the solid. The presence or absence of particle diffusion does not sharply distinguish solids 
from liquids of point particles for similar reasons. Note also that labelling of vortices 
becomes ambiguous when the lines pass within a coherence length of each other — one 
must sum over "direct" and "exchange" connection possibilities to deflne the statistical 
mechanics. Entanglement also has a quantitiative meaning for dynamics: if lines can only 
cross or recombine by overcoming large free energy barriers, such lines are dynamically 
entangled, similar to a dislocation tangle in a work hardened metal. In this paper we will 
continue the usage in the recent literature of referring only to phases which are "superfluid" 
in the sense of Eq. (1.6) as entangled. 

Figure 4 illustrates two distinct scenarios for phase diagrams of vortices with increasing 
temperature, which we call "Type I" and "Type 11" melting. In Type I melting a flrst- 
order transition separates a line crystal with 7^ from a flux liquid with i/jq ^ which 
may or may not be hexatic. In Type II melting, both order parameters are nonzero in an 
intermediate "supersolid" phase. As discussed in Sec. 4, vacancies or interstitials enter 
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the Abrikosov flux lattice at in much the same way as the flux hues penetrate the 
Meissner phase at Hci. Although this transition can in principle be continuous, even in 
the presence of strong thermal fluctuations [17], the melting of the supersolid into a liquid 
is likely to remain a first-order transition, since, as far as is known, continuous melting 
transitions in which materials lose crystallinity in two or more directions do not occur in 
three dimensions. 

Can a real supersolid phase exist in quantum crystals? According to the review by 
Andreev [23], "the experimental data available at present show that the possibility ...(of a 
supersolid)... can hardly take place in [bulk] ^He crystals." In thin films of ^He, supersolid 
phases can, however, exist. On substrates in which one or two layers of incommensurate 
solid ^He form, there can be a regime in which the atoms in the next partially filled layer 
are superfluid; this is a 2D supersolid phase even though there are almost no vacancies or 
interstitials in the close-packed solid layers. In 2D electron crystals at zero temperature, 
some kind of solid phase in which interstitials proliferate may occur [26b] . However, since 
electrons are fermions, such a phase would probably not be supersolid, expect perhaps at 
extremely low temperatures. 

In the case of flux lines, the softer (longitudinal) nature of the interaction between 
lines makes a supersolid less unlikely than in bulk ^He. However, in the limit of continuous 
flux lines, where the boson analogy is best, we shall see that such a phase is still improb- 
able. The new ingredient in the flux system, which does not have an analog for bosons, is 
the discreteness of the layers, i.e. discreteness in the time-like direction. As we shall see, 
this will definintely cause a supersolid for sufficiently large B in strongly anisotropic lay- 
ered superconductors. The characteristic magnetic field above which the layering becomes 
important is the crossover field [12, 24] 

B.-l'^. (1.6) 

the effective mass anisotropy. Physically, this field represents the crossover between differ- 
ent behaviors for the energy at a "jog" where one line is shifted by an intervortex spacing 
^0 = (0o/-B)^/^ from one layer to the next. A jog with shift distance a < do/^ costs an 
energy quadratic in a, while for a > do/j, the energy cost will be linear in a and a jog 



where do is the layer spacing and 
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will tend to spread out over more than one layer. For jogs with shifts oq, the former will 
obtain for B » B^- In a sufficiently anisotropic system, B^ <S Hc2{T = 0) and a regime 
exists, B ^ Bx , in which the vortex fluctuations are predominantly two dimensional. 

Figure 5 summarizes our conclusions about the phase diagram as function of mag- 
netic field B (i.e., vortex density) and temperature T, with external field H is parallel to 
the c axis. In a sufficiently anisotropic system, the supersolid must exist for B ^ B^ , 
where vortices in different copper-oxide layers are approximately decoupled by thermal 
fluctuations [12,24]. The defects in this regime are equivalent to the "unbound disloca- 
tion quartets" discussed by Feigel'man et al. [13]. The existence of a supersolid for fields 
B ^ B X is a more delicate matter, which we shall discuss more quantitatively in Sec. 4. 
For fields B ^ Bd = (po/X^, so that vortices interact via a logarithmic potential, the 
preferred defects arc intcrstitials instead of vacancies. Similar conclusions were reached by 
Fisher et al. and by Cockayne and Elser for electrons interacting with a 1/r potential in 
two dimensions [25]. At low fields {B ~ -Bci), we expect that vacancies are the favored 
defect, as in most crystals with short-range interactions. Although a supersolid is possible 
in principle for any field B < Bx , the numerical estimates of Sec. 4 indicate that B ~ Bx 
is required for this new phase to exist in high-Tc superconductors. 

Even if an equilibrium supersolid does not exist for B < Bx , vacancies and interstitials 
may still appear for nonequilibrium reasons. If, for example, vortices pass through a flrst- 
order freezing transition upon cooling under conditions of constant B, the system will 
initially phase separate into crystallites coexisting with a liquid of a different density. 
When the liquid disappears at a lower temperature, the crystal density must change to 
keep the macroscopic B fleld constant. In the presence of pinning, kinetic constraints may 
force the crystal to absorb interstitials or vacancies rather than to changing its overall 
lattice constant. If freezing occurs above the fleld corresponding to maximum melting 
temperature in Fig. 5, a nonequilibrium concentration of interstitials would then result 
[26]. 

In Section 2 we derive simple estimates for the defect unbinding temperature in various 
regimes, and compare these to the melting temperature. Numerical calculations of the 
energies of various kinds of straight defect lines in the high fleld regime are presented in 
Section 3. In Section 4 these energies are used to estimate Td more quantitatively for Bci < 
< B ~ Bx- We also discuss the consequences of supersolid order for neutron diffraction 
and theories of the melting transition. And, flnally the implications of a supersolid phase 
for resistivity measurements, neutron diffraction and dislocation conflgurations. 
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2. ESTIMATES OF DEFECT UNBINDING TRANSITIONS 



2.1 Model Parameters and Field Regimes 

To obtain the important physical parameters in a thermally excited vortex lattice, it 
is instructive to examine one representative fluxon in the confining potential provided by 
its surrounding vortices in a triangular lattice. 



dz ) . (2.1) 



As we shall see, this simple model gives predictions for melting equivalent to those ob- 
tained via "nonlocal" elastic constants and the Lindemann criterion [27]. Closely related 
models (such as Eq. (1.2)) will be used to determine when defects proliferate. The effective 
potential Vi[r(z)] in (2.1) arises from a microscopic pairwise interaction V^[ry(z)] between 
two parallel flux lines ri{z) and rj{z) with separation rij{z). In the London approximation, 
this potential is [1] 

Vinj) = 2soKo{rij/Xi_) (2.2) 

where 

£o = (</>o/47rAi)2 , (2.3) 
with A_L the in-plane London penetration depth is the energy scale per unit length. The 

1 /2 

Bessel function Ko{x) ^ ln{x) for a; ^ 1, and Ko{x) ^ (^) for x large. The 

parameter si is the tilt energy of a flux line, given approximately by 

si ^ 7^60 H^o/U) {B > Sci) (2.4a) 

where 7^ = M^/M^ ^ 1 is the mass anisotropy, ao = {(j)^/ B)^/"^ , and is the in-plane 
coherence length. Equation (2.4a) follows from the wave-vector-dependent tilt modulus 
evaluated in the short distance regime relevant for melting [24]. This result applies only 
when B ^ Bd = (^o/-^j_7 so that we may neglect the electromagnetic coupling between 
Cu02 planes. When B ~ Sd, this coupling is important and the tilt energy becomes 
[12] 

ei ~ eo , {B ^ Bel) . (2.4b) 

Although the tilt modulus has more generally a complicated wavevector dependence 
[28], it is approximately constant over the (short distance) length scales of interest to 



8 



us here. "Nonlocal" (in z) contributions to the interaction potential [29] are similarly 
unimportant provided [30] 



d-r 



dz 



) <7" 



(2.5) 



which is well satisfied throughout the crystalline phase. As discussed by Brandt [31], the 
same criterion justifies the neglect of higher order terms in [df{z)/dz]'^ in Eq. (2.1). 

Three important field regimes for fluctuations in vortex crystals are easily extracted 
with this approach. We first rewrite this imaginary time path integral as a quantum 
mechanical matrix element [32], 



Z(f^,0;L) = (fl|e-^^/^|0) 



(2.6) 



where \0) is an initial state localized at 0, {f±\ is a final state localized at r±, and the 
"Hamiltonian" 7i is 



n 



2£i 



^ vi + Vi(fi) 



(2.7) 



The probability of finding the fiux line at transverse position fj_ within the crystal is 
i^oirj.), where 0o(n_) is the normalized ground state eigenfunction of (2.7) [33]. 



When B » Bci, the potential is logarithmic, and we expand Vi{f±) about its mini- 
mum at rj_ = to find 

r t2 „ 1 „i 

•00 = SQtpo (2.8a) 



rp2 1 



where (neglecting constants of order unity) 



k 



an 



r=ao 



(2.8b) 



and ao is the mean vortex spacing. Equation (2.8a) is the Schroedinger equation for a 
two-dimensional quantum oscillator, with fi ^ T and mass m — > £i. The ground state 
wave function is 



V'o(r_L) 



1 



with spatial extent 



27r 
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(2.9) 



(2.10) 



Melting occurs when r* = c^ao, where cl is the Lindemann constant, so the melting 
temperature is 



Tm = cl^fl^x ao , (Bel < 5 ^ Sx ) (2.11) 

in agreement with other estimates [27]. Note that Eq and ii are themselves temperature 
dependent so that Eq. (2.11) is an implicit expression for Tm{B), or, equivalently, a melting 
field Bm{T). 

This result is only applicable when ^ Hc2{Tm) = (/>o/[27r^^(T^)] where C±{T) 
is the temperature dependent in-plane coherence length including the efi^ects of critical 
thermal fiuctuations which make Hc2(T) ~ (Tc — T)^/"^ for T sufficiently close to the 
zero-field transition at Tc. In this regime A_l ~ (Tc — T)~^/^ and the resulting ratio 

Bm/Hc2{Tm) = const. 

independent of magnetic field or anisotropy. The logarithmic factor in £i, Eq. (2.4a), is 
thus of order unity since ao/^±{Tm) = const. Such critical effects were not taken into 
account in calculations of the melting temperature prior to Ref. 12. 

Vortices in the crystalline phase will travel a perpendicular distance r]_ in a "time" 
along the z axis where [7] 

^ rl/{T/e,) 

% «o • (2.12) 

A new high field regime arises when ~ where do is the average spacing of the 
copper-oxide planes, i.e., for B ^ B^, with [12,24] 

Bx ~ Ji ^ , (2.13) 

So Gq 

which is up to a logarithmic factor the same criteria discussed in the previous section, 
Bx ~ (7^0o/<^o)- Well above this field, the planes are approximately decoupled, and 
Tm may be estimated from the theory of two-dimensional dislocation mediated melting 
[34,12,24] 

Tm ~ 0-5|^ , {B:>Bx). (2.14) 
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where the prefactor 0.5 is a rough estimate of the effects of phonon fluctuations [34]. Note 
that the convergence to Eq. (2.14) for high fields may be quite slow due to the very strongly 
divergent translational correlation length at the melting temperature in two dimensions 
[12,24]; this causes T^(-B) to appraoch Eq. (2.44) only as l/ln^(B/5x). 

The estimate (2.11) also breaks down at low fields B ~ Bd where the logarithmic 
interaction potential must be replaced by an exponential repulsion. The two-dimensional 
harmonic oscillator approximation can again be used with the replacement 

k^^e-''^/^^. (2.15) 
The transverse wandering distance is now 



2x2 \ 1/4 



and takes place over a longitudinal distance 



i \ ^ao/4A_, 



(2.16) 




A±e""/2^^. (2.17) 



The low field melting temperature becomes 

Tm ^ civ^ ^e-"°/^^- {B ^ Bel) , (2.18) 

consistent with earlier predictions [7]. Although we have retained the distinction between 
£i and Eq in these formulas, note from Eq. (2.4b) that Si ~ Sq in this regime. 

The predictions (2.11), (2.14) and (2.18) are combined to give the reentrant phase 
diagram for melting shown in Fig. 5 [35,12]. We note that when fluctuations are rela- 
tively weak, so that melting occurs closer to Tc, Eqs. (2.11) and (2.18) must be solved 
self-consistently for using the temperature dependence of cq, ei and Aj^. This proce- 
dure causes the maximum in the melting temperature to bend downwards towards in 
materials like YBCO [12]. B rather than H is used for the vertical axis since it is B rather 
than H which is flxed for the thin, low aspect ratio crystals usually studied in experiments. 
This flgure is only schematic because a flrst-order transition [2] at constant B would lead to 
two-phase coexistence with domain size limited by the strong interactions between vortex 
tips as they exit the sample surface in thin crystals. Analytic estimates and boundaries 
for melting in the various regimes are summarized in Table 1. 
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2.2 Vacancy and Interstitial Unbinding Transitions 

Consider a defect such as the vacancy shown in Fig. 2. To estimate when its free energy 
changes sign, we need to determine the kink separation which appears in Eq. (1.3b). 
Up to factors of order unity, a similar estimate for would apply to an interstitial vortex 
which, for example, hops between the centers of the triangular cells of an Abrikosov crystal. 
Interstitial wandering would then take place on a honeycomb lattice, and the factor In 6 
in Eq. (1.4b) would be replaced by In 3. 

We suppose initially that B^x <^ -B ~ 5x, so that the pair potential is logarith- 
mic with well-coupled Cu02 planes. As shown in Sec. 3, the defect line energy is then 
proportional to the characteristic line energy scale £o of the pair potential (2.2), 

ed = CdSo . (2.19) 

In Section 3, we determine the constant Cd for various kinds of defects, and show that in 
fact interstitials, rather than vacancies are favored energetically for B ^ Bd- 

An interstitial surrounded by a triangular cage of nearest neighbors will fluctuate in 
much the same way as the vortex in the hexagonal cage discussed in the previous section. 
To "tunnel" to one of its three neighboring triangular plaquettes, it must overcome a 
barrier of order sq- The energy associated with one of these kinks or tunneling events is 
of order 

Ek ~ ao , (2.20) 

from Eqs. (2.12) and (2.1) since the "kinks" are spread over a distance This expression 
is just the WKB tunneling exponent which follows from the "quantum mechanical" analogy 
discussed above. The line-density of kinks is then of order 

nkink ~ ^ e-^'=/^, (2.21) 

where £^ should be the same order of magnitude as in Eq. (2.12). Similar estimates 
would apply to vacancies: Kinks arise in Fig. 2, for example when one of the flux lines 
surrounding an unoccupied site jumps over a barrier in the lattice potential Ui{r) of order 
Eq and essentially changes places with the vacancy. The "attempt frequency" associated 
with this event should again be of order 
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For both types of defects, we conclude that the spacing between kinks is of order 
= K^^, or 

4 ~ ^°e-^'=/^ . (2.22) 

Upon substituting Eqs. (2.19) and (2.22) into (1.3b), we find a simple self-consistent 
equation for T^, 

Td = ciy/eoSi ao exp[c2 V^o^iOo/Td] (2.23) 
where ci and C2 are constants of order unity. It follows that 

Td = C3^/^ao (2.24) 

which has the same form as Eq. (2.11) but with C3 a different numerical constant. Therefore 
a supersolid phase will occur whenever C3 < c|. We attempt to estimate C3 in Section 4. 

When B » Bx, defects can hop one lattice constant or more between copper-oxide 
planes of spacing do, and the defect path integral (1.3) must be evaluated in a different 
limit. We now neglect the lattice potential Ue{fd) entirely and evaluate the functional 
integral which remains by discretizing the path integral along z in units of do. The result 
is 

which leads to a defect unbinding temperature 

Td{B) ^ const. X eodo/ In[27r5/Sx] . (2.26) 

Since Td ~ l/ln(27rS/S x) in this regime, while is asymptotically S-indcpendent 
according to Eq. (2.14), defects must proliferate for B ^ Bx- Equation (2.26) agrees with 
the unbinding transition of dislocation pairs from bound dislocation quartets estimated in 
Ref. 13. 

The result Eq. (2.26) for the limit Bx ^ B can also be obtained from a simple physical 
argument. In the limit of nearly decoupled layers with a very weak Josephson coupling, 
the density of defects in each layer is just Ud ~ e"'^'*/^ where 

Ed = doSd , (2.27) 

is the energy of the defect in a single layer. When the interactions between layers are weak, 
they can be estimated perturbatively. If a vortex passing through two neighboring layers 
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is misaligned by an amount 6, the excess cost is of order ^'^sob/do from the Josephson 
couphng between the layers. Naively, the cost of misalignment of the vortex which makes 

— 1/2 

up an interstitial "line" is obtained by setting b = aon^ ' , the typical jog distance of the 
defect between nearly decoupled layers. This is certainly an upper bound for the kink-free 
energy since fluctuations of the other vortices in each layer will add a large incoherent part 
to the phase difference associated with the interstitial. We conjecture that this will result 
in a reduction of the full Josephson energy to 

E, ^ (2.28) 

"0 

with A; < 1. The free energy per layer of the defects is then 

Fd ^ EdUd + Tud In nd + ndEi{nd) . (2.29) 

Minimizing with respect to yields the result that Fd becomes negative, and hence Ud 
positive, when T > Td with 

Td^kEd/HB/B^) (2.30) 

just the form of Eq. (2.26). Note that the magnetic interlayer coupling will only add an Ud 
independent term to Ej because of the screening due to relaxation at the other vortices as 
discussed in Appendix E. 

The 2D defect energy Ed is correctly obtained via Eq. (2.27) from the 3D calculations 
of Ed described in this paper since in the absence at Josephson coupling the vortices in 
2D layers interact logarithmically at all distances with the effective penetration length at 
long distances only reduced by negligible O{do/X±) corrections for small do. 

Finally, when B ~ Bci , the exponential interaction between vortices leads to a differ- 
ent result. The lowest energy defect in this regime of short-range interactions is presumably 
a vacancy, as is usual for solids with short-range interactions. The characteristic defect 
line energy is now 

Ed = const. X -f £oe"''°/^^ . (2.31) 
^± 

The barrier to produce a kink in the trajectory of a vacancy or interstitial has similar 
exponential dependence on ao, so the resulting kink energy is 

Ek ~ V^i ^ e-«°/2^- . (2.32) 
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Upon using Eq. (2.17) for £° in Eq. (2.22), we find via Eq. (1.30) that 



r T 

Fd = SdL 1 — const. X — e 



Ek/T 



(2.33) 



which (provided the constant is not too small) changes sign for 



Td ~ const. 



const. X \^7oii 




0_ ] -ao/2X^ 



(2.34) 



in the limit of low fields. A supersolid phase is in principle possible for low fields if Td < Tm 
as given by Eq. (2.18). We conclude in Sec. 4 however, that < Td for B ^ B^, which 
makes this possibility rather unlikely. 

Our estimates for interstitial or vacancy proliferation temperatures are summarized 
in Table 1. Combining these estimates for Td and (and assuming that the supersolid 
is unfavorable for small fields) leads to the schematic phase diagram shown in Fig. 5. 
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3. FORMATION ENERGY OF LINE DEFECTS 



In this section we discuss more precisely the defect line energy for rigid configu- 
rations of straight vortex lines. Numerical calculations are presented for various types of 
defects in the limit Bci B <^ Hc2- 

3.1 Gibbs Free Energy of a Perfect Vortex Lattice 

Let x{t) be the equilibrium position of the £th flux line. In the London limit, k = 
^±/^± ^ 1) the pair potential between two straight fliix lines is given by Eq. (2.2), 

VM = ^§yTKo{rw/\A.) , (3.1) 

where rw = \x{£) — x{£')\, and Kq is the modified Bessel function, with the asymptotic 
behavior Ko{x) ^ (Tr/2x)^/'^e~^ for large x and Ko(x) — ln(a;/2) — 7 for small x. 

For fl^ II S II c the Gibbs free energy per unit length L and unit area A is 

g = n(^ei- ^^^^ + ""^oXl ^o{ro£/>^±) ■ (3.2) 

Here ei = (^j^xi) In/? = eo Inre is the energy per unit length of a single flux line (ignoring 
an additive constant correction to Iuk), and n = N/A is the number of lines per unit area. 
The second term in the Gibbs free energy represents the interaction energies of the lines, 
while the first includes the effect of the external magnetic field H and favors large values 
of B. H plays the role of a pressure (or chemical potential), which tends to increase the 
density of lines. The prime denotes that self-interactions (the £ = term) are omitted 
from the summation. The equilibrium magnetic flux density i?eq = ncjiQ is obtained by 
minimizing the Gibbs free energy, i.e., by solving dG/dB = 0. In the limit i?eq S> 4>o/\']_ 
one flnds [36] 



+ 0{B;') , (3.3) 



where Hd = The constant C = — ln47r — 2 + 7 + ^ depends on the lattice structure 

with 7 fa 0.5772.. is Euler's constant. For a hexagonal and square lattice we flnd Aq fa 
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0.0797107, and A4 0.1008797, respectively [36]. For a rectangular lattice with lattice 
constants Lx and L„ the value of A has to be calculated from 



R 



LxLy 



+E 



exp 



-CLxLy/ATT 



G^LxLy/AiT 



(3.4) 



The corresponding minimum value of the Gibbs free energy is 



(3.5) 



We arc interested in the formation energy of "point defects" at constant line density n 
corresponding to a fixed external magnetic field H. This later condition allows us to 
consider only the interaction part of the free energy (3.2) for N particles, i.e.. 



En = Neo^Koiroe/X±] 



(3.6) 



in determining defect energies. With Eqs. (3.2)-(3.5) one finds 



En = 2Neo { nnXl + ^ 



In 



1 



nXl 



+ C + 1 



(3.7) 



3.2 Definition of Defect Energies 

3.2.1 Constant lattice spacing 

We follow Ref. 25 and illustrate the definition of defect energies for the case of a 
vacancy line. Consider first a perfect lattice with + 1 fiux lines confined to an area A. 
The corresponding interaction energy is En+i- Removing a line at the origin without 
allowing the flux lattice to relax reduces the interaction energy to 

E' = ^E'^(^o£) - E'^(^o.) ' (3-8) 
e e 

i.e., we have subtracted the interaction energy of the removed center line with the rest of 
the lines. Hence the (unrelaxed) defect energy for a vacancy at constant lattice spacing is 
given by 

^ ^/ _ ^^^^ ^ __1_E^^, . (3.9) 
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3.2.2 Constant line density 



The vacancies that occur in a superconductor occur at constant hne density or constant 
chemical potential, not constant lattice constant. Starting with a large perfect lattice 
containing N flux lines in a fixed area A one can imagine rearranging the flux lines (with 
N fixed) in such a way that the resulting configuration is identical to the one generated 
above by removing one line from a perfect lattice containing N + 1 lines. Thus, the 
unrelaxed vacancy energy at constant line density is defined by [25] 

E^ = E'-En = + (En+i - En) . (3.10) 

Upon using Eqs. (3.2)-(3.7) one gets for the Gibbs free energy per unit length of a vacancy 
line in an unrelaxed lattice. 



Y = Ey- — 



In ( ]+C + 2 



(3.11) 



Note that we have neglected terms of the order n/N = 1/A since we are interested in the 
thermodynamic limit with n fixed and N oo. The defect energy itself is ey = Gy ~ 
Co In + const, which is the correct order of magnitude except for the logarithmic 

divergence as oq/Ax —>■ 0. As we shall see, this singularity disappears when we allow the 
lattice to relax around the defect. 

3.3 Elasticity Theory 

Before allowing static relaxations, we review the energies of the phonon distortions 
which are involved. When each vortex line is given an arbitrary displacement ■u{£) from its 
equilibrium position x{£) the interaction (potential) energy of the vortex crystal is given 

by 

^=^E'^(I^(^^')+^(^^')I), (3.12) 

where x{££') = x{£) —x{£') and u{££') = u{£) — u{£') are the differences in the lattice vectors 
and displacement vectors, respectively. Expanding up to second order in the displacement 
fields yields 

£; = E^ + -L J]J]0«/'(£,f)z.-(£)t./3(f)+ (3.13) 



2A 

i,a. £',13 
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where Ej^ is the interaction energy of a rigid flux lattice. Here a,j3 = x,y label the 
Cartesian components and the dynamic matrix is given by 



Upon defining the Fourier transform 



32 



for e^e' , 



with 



ttr lim 



x^o dx"dxf^ 



the interaction energy can be written as 



(3.14) 



C"'3(g) = ac 5^ (/)"^(£, 0)e-^^"'-^(^'°) = - (5"^(g) - ^"^(0)) , (3.15) 



, (3.16) 



(3.17) 



Here Qc = = n ^ is the volume of the unit cell of the triangular lattice and we have 

introduced the short hand notation J^ — J integrals over the in-plane wave vector 



In the limit of a very large London penetration depth, Aj^ 3> Oq, the dynamic matrix 
takes the form (see Appendix A) 



(3.18) 



for small values of q. Upon comparing with the usual continuum elastic description of 
vortex sohds [27], C"^{q) = [cil{q)q°'q^ + CQ6{q)d°'^q'^], we find that the bulk and 
shear moduli are 



cii(g) fa n 



2 47reo neo 



(3.19) 
(3.20) 
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where n = is the vortex density. In this hmit of rigid, parallel vortices we obtain no 
information on the ^fz-dependence of these quantities or on the tilt modulus C44. Equa- 
tion (3.20) agrees with the shear modulus first obtained by Fetter et al. [38]. Equa- 
tion (3.19) is correct for A^^ < qj_ < clq^, which encompasses a wide range of length scales 
when B > (po/X^. 

3.4 Variational Calculation of the Relaxational Energy of a Straight Vacancy Line 

In this section we calculate the relaxation energy of a straight vacancy line by a 
variational approach. The presence of a vacancy causes a distortion of the lattice described 
by the displacement field {!{£). Assuming the displacements ■u{£) to be small and slowly 
varying the relaxational energy for a vacancy can be obtained by minimization with respect 
to u(e) of 



£1' 



-1 J] Vl« (£)«"(£) 

-lj2V2^£>''i£>^i£), (3.21) 



where 

Vi°W =-isV(x) (3.22) 

and the dynamical matrix (f)'^^ {£,£') is defined in Eq. (3.14). 

The first term in Eyk{{u}) is the elastic energy of the lattice distortion caused by the 
vacancy. It is overcounts the energy of the missing center line with the rest of the crystal. 
This contribution is subtracted by the second and third term in Eq. (3.21), which is just 
a Taylor expansion of J^e^^o ^(1 ^(■^) + '^(•^) ~ ^(0)1)' where ^(0) = is the position of the 
vacancy in the lattice. 

In Fourier space, we have 

' (g) {£)e~'^'^^^^ (3.24) 



u 
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C-^i^ =a, J2 Mi, O)e-^^^--^(^) (3.25) 

i 

=a, J2 V^^{e)e-'^-^^^'^ (3.26) 

£=^0 

V^^'^iq) =ac J2 {tfe-''^-^^^^ (3.27) 
and the vacancy relaxation energy becomes 

J a 

y2^{?- k)u''{k)u^{-q) . (3.28) 



1 




2ac Jq Jk 

Following Ref. 25, we make a variational Ansatz for the lattice distortion field, 

u^i^ = iac^fiq), (3.29) 

where f{q) = 1+cq+dq^ with c and d as variational parameters. As discussed in Appendix 
B, the constraint /(O) = 1 is enforced by the long range potential. Ff(^ and V^^{q) can 
be expressed in terms of the dynamic matrix C°'^{q} 

Vrio) = ac E Vr{e)e-^'^<'^ = -- i^Vr{q) , (3.30) 

F2«^(^ = acJ2 V^^tie-'^-^^^^ = -C'f^iq) + V^'''^{q = 0) , (3.31) 

where we have used the asymptotic form of the potential, limjj^j^^co y(a;) = —\nx, to 
derive equation (3.30). This gives 

¥{-{^ = 2180 (3.32) 



Vg2 2 

and for the above Ansatz for the displacement field one can take 

V,"^{q) = -C-^iq) . (3.33) 

Upon using the long wavelength approximation (2.18) for C°'^{q), one finds 

E^iiM^)}) = 2^0 {h + I2+ h) (3.34) 
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with /j given in Appendix B. The variational calculation for c and d yields (see Appendix B) 

1 



d = 



7kd 
_ 10 



0.161^0;, 



-0.114ac , 



(3.35) 



(3.36) 



where kd = •\/47r/ac is a cutoff which preserves the area of the underlying hexagonal 
Brillouin zone. Hence, the elastic relaxation contribution to the vacancy energy density at 
constant lattice spacing is 



^vR-y 



In 



1 



+ 



265 

252 



(3.37) 



The total vacancy energy density at constant density instead of constant lattice con- 
stant, is, following the considerations of Sec. 3.2, 



G 



V 



(3.38) 



The last three terms are given by Eq. (3.11), which, when combined with (3.37), leads to 



G 



V 



-7 - ^6 + 



265 
252 



(3.39) 



Our variational vacancy line-energy is thus ey = Gy, or 



ev = 0.1973£o 



(3.40) 



Note that the logarithmic divergences as ao/X± — > in Eqs. (3.11) and (3.37) have cancelled 
to yield a finite result. 

A heuristic "back of the envelope" argument for the vacancy energy can also be con- 
structed: assume that the phonon displacement in real space for a vacancy at the origin 
takes the isotropic form 

(3.41) 



27r x^ii) ' 

consistent with a six-fold symmetry of this defect site [37]. The parameter ^Iq is the 
area change induced in the flux line lattice. To keep the density of flux lines fixed (as is 
appropriate for B ^ (f)o/X'j_), we take Qq = ac, the area of one unit cell to cancel the 
vacancy energy. Since V • tt = and dgU = 0, the only contribution from the continuum 
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elastic free energy per unit length comes from the wave-vector independent shear modulus 
term, 

= J d^xceeu^j , (3.42) 

where Uij = | {diUj + djUi), i = 1,2 and [38] cqq = | no£o- In Fourier space, this expression 
becomes ^ 

ey = C66 j ■^^\uij{q)\'^ , (3.43) 

with Uij{q) — —ttcQigj /q'^ . Upon using Eq. (3.20) and imposing a circular Brillouin zone 
of radius ka = ^/ATr/ac, we find finally 

ev = £o/4 , (3.44) 

an estimate only 25% greater than Eq. (3.40), and very close to the numerical value for 
the relaxed vacancy configuration with six- fold symmetry [37] . 

3.5 Numerical Calculation of the Defect Energies 

Interstitial defects generally occupy lattice sites of low symmetry, so that analytic 
methods become quite laborious. In this section we describe numerical calculations of the 
various defect energies which do not require the approximations used above. Our goal is 
to calculate the defect energies for an infinite system. 

In the computer simulations, we have not an infinite system but a system with a large 
but finite block of particles with periodic boundary conditions. In order to handle the long- 
range logarithmic interaction between the flux lines we use the the Ewald sum technique, 
which amounts to including the interaction of the flux line with all its periodic images [39]. 
For comparison with the computer results, one must therefore calculate not the energy 
of a single defect in an infinite system, but the energy per block of an periodic array of 
defects in an infinite system. After correcting for these "image" defects one can extract 
the desired energy of a single defect. In practice we choose the system large enough that 
the interaction of the defects with its periodic images can either be neglected or calculated 
by means of linear elastic theory. 

The interaction energy of the fiiix line lattice per unit length in the z direction is 

£^=^E^o(r«VA±), (3.45) 
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where our energy unit is £'0 = 2£o- This sum over the infinite lattice can also be written 
as 

E <^(^^^'/A±), (3.46) 

1,1' e box 

where the summation is over all particles within the basic box of size (La;, Ly). The effective 
pair potential within the basic box 

<i>{rii'/\i_) = E' ^o(| m') + R I (3.47) 

R 

represents the interaction energy between a flux line at position x{t) and one at position 
x{l') together with all its images at positons x{l') + R. The sum over R runs over all 
simple cubic lattice points, R — {£Lx,mLy) with £, m integers. This vector reflects the 
shape of the basic box. The prime on the summation sign of Eq. (3.47) indicates that we 
omit the term ^ = for x{U') = 0. 

For a numerical calculation equation (3.47) is not a suitable starting point because of 
its poor convergence properties. Therefore, we make use of Ewald's summation method 
[39]. Using the integral representation of the modified Bessel function Kq and Ewald's 
generalized theta function transform (see Appendix A) we find for ^ in the limit of 
large A 




e 



iG 



5.- exp [-L^LyG'^6/4TT] 



2^ L^LyG^/AiT 

G 

where Ei{x) = dyy~^e~y is the exponential integral function. The vectors G = 
27r{m/Lx, n/ Ly) with m and n integers index the square lattice reciprocal to the lattice of 
image lines. Note that this result is valid for any choice of the Ewald separation parame- 
ter 5. It can be used as a numerical check and to optimize the convergence properties of 
Eq. (3.48). We choose 5=1. 
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For 0(0) we get 



R 



+ 



7 



L^Ly 2 ' 2 47rA^ ^ 2 



27rA2 1 , 
^ + - In 

LxLy 2 



^)+^(7-l-ln47r + ^rect) 



(3.49) 



The value of the Ewald sum A^-ect depends on the shape of the basic box (compare with 
Eq. (3.4)). We choose the rectangular basic box such that L^^ = 5ao • t and Ly = 3v^ao ■ t 
with t integer. Then one gets from a numerical evaluation of the Ewald sums Arect ~ 
0.1018412 and the interaction energy in units of Eq becomes 



E= ^ 

2 LxLy 



+ 



N 



7-l-ln47r + ^rect + ln 



LxLy 



+ \Y^^tu'/\a.) (3.50) 



with the effective pair potential 



1 



{x-Rf-KX 1 



LxLy 



-iG-x 



exp [-LxLyG'^/iir] 1 
LxLyC^/AiT 2 



(3.51) 



R ^ ' G 

Upon comparing Eq. (3.50) with the interaction energy of a perfect fliix line lattice, 
Eq. (3.7), one gets for a rectangular basic box 



1 / ^ jv 



(3.52) 



where A?" = 30 • is the number of flirx lines in the box of size [Lx-, Ly) = (5, 3-\/3)too. 

For the numerical calculations it is sufficient to consider only those parts of the in- 
teraction energy, which explicitely depend on the coordinates of the flux lines. Hence we 
consider the quantity 



£,£' 

The force between flux lines at a distance x{££') is given by 

Fixer) = Fee = 'g^^^M ■ 

Using the above expression for the interaction potential $ one can write 

exp — \ X — R\^ tt/LxL^ 



(3.53) 



(3.54) 



^/ 1 \ ^ 2tv 



I X — -R p n/LxLy 
^^'n ■ (n -\ exp {G'^LxLyjAm] \ 



(x-R>^ 



(3.55) 
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These forces are most efficiently calculated together with the energy (in one subroutine). 

We now explain the algorithm by which the formation energies of the various types of 
defect lines are determined. We start with an initial configuration {x^^\i)} which, after 
relaxation, will contain the desired defect. This is easily achieved by adding or removing 
lines from the perfect hexagonal lattice. This procedure leads to the defect energies at 
constant lattice constant, and hence must be corrected as in Eq. (3.38) to produce energies 
at constant vortex density. The corresponding initial configuration for a vacancy and a 
variety of other defects are shown in Fig. 6. 

The lattice relaxation process was performed by standard methods adapted from 
molecular dynamics simulations [40]. For advancing the positions and velocities of the 
flux lines we implemented an artiflcial dynamics via a leapfrog algorithm 

Mt + l^t) = Mt - l^t) + stMt) (3.56) 

Xiit + St) = Xiit) + St nit + ^St) (3.57) 

The acceleration a£ — fi/moi the £th flux line is calculated from the forces = 'Y^ii^i F^i. 
In each step of the iteration we calculate the kinetic and potential energy, and the forces 
using Eqs. (3.50), (3.51) and (3.55). This procedure is repeated until an equilibrium 
configuration is reached. The "mass" m and the time step St were chosen to accelerate 
the convergence to equilibrium. Equilibrium here means that the forces on the flux lines 
become zero. Hence the method is capable of flnding not only minima but also some saddle 
point conflgurations, at least for initial states with high symmetry. 

In order to test the accuracy of our method we calculated the energy of a perfect fliix 
line lattice with AT = 30 x particles for t = 1, 2, 7, and compared them with the exact 
result Eq. (3.52). Our results are summarized in Table 2. The relative difference is less 
than 10~^ and can mainly be attributed to the inaccuracy in the numerical approximation 
we have used for the exponential integral function [41] . 

Starting from the initial conflgurations shown in Fig. 6 we have determined the relaxed 
conflgurations using the algorithm described by Eqs. (3.56) and (3.57) for various system 
sizes. From this analysis we can extrapolate to the defect formation energies of an inflnite 
system. 
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For the vacancy we find that starting from the initial configuration in Fig. 6, which 
has the six-fold symmetry of the hexagonal lattice, the lattice first relaxes into a saddle 
point configuration, which possesses the full symmetry of the lattice. This configuration, 
however, is unstable with respect to a compression along one of the three axis connect- 
ing the nearest neighbors at the vacancy. It finally relaxes into a configuration of lower 
symmetry (see Fig. 7). Note, that there are three equivalent orientations of this relaxed 
vacancy configuration. In Table 3 we summarize the formation energies of the stable and 
saddle point configuration of the vacancy for different system sizes. 

The edge interstitial relaxes starting from the initial configuration in Fig. 6 into a 
saddle point configuration shown in Fig. 8. From the numerical simulation we find that 
this configuration is unstable with respect to small amplitude "buckling" perpendicular to 
the edge. The edge interstitial relaxes into a "buckled configuration", which is identical 
to the relaxed centered interstitial configuration, also shown in Fig. 8. These results are 
analogous to findings by Cockayne and Elser for the two-dimensional Wigner crystal [25], 
where the edge interstitial is also unstable with respect to "buckling" perpendicular to 
the edge of the triangle. However, the time required for this relaxation process is much 
larger than the relaxation time from an initial edge interstitial configuration constrained 
by symmetry to go to a final edge state. Thus, the interstitial appears to occupy a very 
flat minimum in configuration space. 

In Table 3 we have listed the defect energies for a vacancy and several types of in- 
terstitials for various system sizes. Whereas there is a clear energy gap between vacancies 
and interstitials, the interstitial energies are all very close. The energy differences are less 
than 1%! The system size dependence for the centered edge, and split centered interstitials 
and two types of vacancies are shown in Figs. 9a and 9b, respectively. 

The lattice conformations resulting from relaxing an edge interstitial and a split edge 
interstitial initial configuration (constrained now not to buckle as in Fig. 8) are shown in 
Fig. 10 (also shown by filled circles is the perfect hexagonal lattice). As can be inferred 
from this figure the configurations essentially differ only by a parallel shift along the edge 
of the triangle. Because the energies are quite close (see Table 3) we conclude that gliding 
of this type of defect along the direction defined by the edge of the triangle (in the absence 
of buckling) must be a low energy excitation. The barrier for this motion is presumably 
of the order of the difference in energy of those configurations, i.e. Agiide ~ 10~^£'o! After 
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very long relaxation times both types of initial conformation, i.e., edge and split edge 
interstitial, may in fact finally relax into the same final configuration. 

We conclude that interstitials, rather than vacancies, are clearly favored at high fields 
and that the centered interstitial is the most likely candidate for producing a supersolid 
in this regime. Since the differences between the energies of the various interstitials are 
so small, an interstitial will have substantial extra entropy, lowering its free energy even 
further. 

3.6 Interaction between Point Defects 

Following Cockayne and Elser [25], we can use defect energies for different system sizes 
to infer the distance dependence of the interaction energies. As explained at the beginning 
of section 3.4, periodic boundary conditions always yield a rectangular superlattice of 
defects. Since we are changing only the size of the big box {L^, Ly) and not its shape, the 
system size dependence of the formation energy should scale the same way as the radial 
dependence of the interaction between two single defects for large and Ly. 

From the system size dependence of the centered and edge interstitial energies, plotted 
in Fig. 9a, one can infer that both defects show an attractive interaction at distances larger 
than 5 lattice spacings. Over the range studied 5 < r < 35 the interaction approximately 
scales as for centered interstitials and as for edge interstitials. From a nonlinear 
least square fit over the limited range we get Ef2t ~ 0.0742 - 0.0787 * N'^-^^ correspond- 
ing to r~^-^^ and Ej*^}; 0.0732 - 0.0078 * N~^-^^ corresponding to r~^-^^ for the edge 
and centered interstitial, respectively. As we shall see, these are definitely not the correct 
asymptotic long distance behaviors. 

We have also analyzed numerically the distance dependence of the interaction between 
two centered interstitials at smaller distances. In order to calculate this energy we have 
started with an initial configuration, that contains two centered interstitials at a distance 
d (measured in units of the lattice constant a) in a perfect hexagonal lattice containing 
N = 480 fiux lines. The interaction energy is found by subtracting from the resulting 
relaxation energy the energy of two isolated single centered interstitials. The results for 
two different directions are displayed in Fig. 11a. If the vector connecting the two centered 
interstitials points along one of the unit vectors of the primitive cell of the hexagonal lat- 
tice, the interaction is attractive up to d ~ 3 a and becomes repulsive for larger distances. 
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In the direction perpendicular to one of the unit cell vectors we find a minimum in the 
interaction potential at a distance of approximately two lattice vectors. (Note that dis- 
tances in Figs. 11 are the distances between the defects in the initial configuration). The 
system size difference discussed earlier indicates that, roughly, some angular average of the 
interaction is attractive for r > 5a. 

In Fig. 9b we have plotted the system size dependence of the symmetric and "crushed" 
vacancy configurations. Whereas for a symmetric vacancy configuration the formation 
energy increases with increasing number of flux lines AT, it decreases for the "crushed" 
vacancy configuration. From a nonlinear least square fit to the data in Fig. 9b we find 
^in? ~ 0.1076 - 0.312 * iV-115 corresponding to r'^-^ and E^Jl ^ 0.125 + 0.0576 * iV-O-64 
corresponding to r~^'^^ over the limited range 5 < r < 35 for the symmetric and "crushed" 
vacancy, respectively. Hence wc conclude that symmetric vacancies appear to have an 
attractive interaction at relatively long distances while "crushed" vacancies, which have 
the lower formation energy repel each other at comparable distances. 

For smaller distances wc have performed calculations analogous to those for the cen- 
tered interstitial. The results for the stable "crushed" vacancy (V2) for the interaction 
along (solid line) and perpendicular (dashed line) to an edge of the unit cell are shown in 
Fig. ll.b. Whereas the interaction energy is attractive for all distances less than 11a along 
the directions perpendicular to the edges of the unit cell, it is attractive for small distances 
and becomes repulsive for distances larger than d ~ 5a for the interaction along the edge 
of the unit cell. This has to be compared with the roughly angular averaged attractive 
interaction for d > 5a obtained from the finite size analysis. 

It is not possible to study the interaction between symmetric vacancies at short dis- 
tances. This is because the anisotropy of the stresses induced by the interactions causes 
the vacancies to deform to the anisotropic crushed vacancy configuration which has lower 
energy. The symmetry axis of the resulting crushed vacancies is parallel to their separtion 
vector. 

It is instructive to compare these results with those obtained from linear elasticity 
theory which should be valid for very large separations. It can be shown that the Inr 
interactions between unrelaxed defects will be completely screened by the relaxation of the 
other vortex lines. This is already evident in the calculations in section 3.4 for a single 
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vacancy: the vortex lines far from the vacancy relax just so as to cancel the overall "charge" 
of the vacancy. The long distance interactions between defects thus have the same form 
as for short-range interactions. These depend on the symmetries of the defects involved. 

For defects with three- or six-fold symmetry — the symmetric vacancies and centered 
interstitials — the interactions are exponentially small in the absence of anisotropy of the 
elastic interactions with the anisotropy associated with the six-fold symmetry [which ap- 
pear at order q'^ in the elastic matrix C°'^{q) of Eq. (3.17)], the interactions will decay as 
cos6^/r^ with 9 the angle between a lattice vector and the inter-defect separation vector. 
The sum over all the inter-defect interactions with the periodic boundary conditions used 
in the numerical computations will thus almost cancel as the systems used are almost 
square with Ly/L"^ = 25/27. The resulting asymptotic L dependence of the vacancy en- 
ergy should thus have a leading term with a small coefRcient whose sign depends on 
details of the q dependence of the elastic matrix and the stresses induced by the vacancy 
which we have not calculated. 

The interactions between defects with only a two-fold symmetry axis-edge and split- 
centered interstitials and crushed vacancies are longer range since they would decay slowly 
even in an isotropic elastic medium. The interactions will have two comparable contribu- 
tions at long distances. 

^ ^if^ + (3.58) 

where is the angle between the symmetry axis of the defect and the separation vector. 
The coefficient K4 is always positive (i.e., repulsive) while the sign of K2 depends on 
the stresses induced by the defects. In the almost square system with periodic boundary 
conditions used in the numerical calculations, the cos 20 term will almost cancel and the 
leading size dependence of the defect energies will be 1 /L^ with a sign cos Aa which depends 
on the cosine of the angle a between the x direction and the symmetry axis of the defects. 
At long distances, however, the actual sign of the interaction between two isolated defects 
will depend on the direction of the defect separation in a different way if 1^2! > K4. 

The apparent distance dependence seen in the numerical calculations is a fit over 
limited range of distance and presumably indicates that the asymptotic behavior discussed 
above has not been reached. 
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4. EQUILIBRIUM DESCRIPTION OF THE SUPERSOLID PHASE 



We now assume that centered interstitials are energetically preferred (as found for 
zero temprature in Sec. 3), and study the low lying excitations via a simple tight- binding 
model. We find a more quantitative estimate of the transition temperature Td^ discuss the 
nature of the transitions at and Tm and compare the properties of the solid, supersolid 
and liquid phases. 

4.1 Tight-binding Model and the Transition Temperature 

As discussed in Sec. 2, the defect proliferation temperature is determined by the 
large L behavior of the partition function (1.2). Equation (1.2) may also be written as an 
integral over a quantum mechanical matrix element analogous to Eq. (2.6), 

Zd = j d^U J d^fj (r/le-^-^^/^lr,) , (4.1) 

where the integral is over entry and exit points fi and fj for the interstitial and Hd 
is an effective quantum Hamiltonian describing the transfer matrix. We adopt a low- 
temperature perspective and use for Tid a tight-binding model which assumes that the 
centered interstitial sits at the sites of a simple honeycomb lattice. As shown in Fig. 13, 
the honeycomb lattice may be viewed as two interpenetrating triangular lattices with sites 
{rA} and {rs}, connected by displacement vectors {Si, i = 1,2, 3} with 

and b = ^ a with ao the triangular lattice constant. 

The interstitial can "tunnel" from one site to another by passing through an edge 
(interstitial) state, which according to Table 3, has an energy nearly degenerate with the 
centered interstitial. The tight-binding Hamiltonian is then 
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Tid = uo 



TA 



I^a) (fA + + ^ Ife) (re + 5b\ 

.rA,5A vbjSb 



(4.3) 



where Ir^) and [re) are normalized states on the two sublattices, {5a} = {Si} and {Sb} = 
{—5i}. Here uq is a site energy, and t is a tunneUing matrix element. 

The eigenvectors of (4.3) are the linear combinations ip±(k) = -^({k, A) ± \k, B)) of 
normalized plane wave states, 



\k,A) 



'No- 

TA 



\KB) 



(4.4) 



where A^o is the number sites in one sublattice and k is confined to a hexagonal Brillouin 
zone. The eigenvalues consist of two bands. 



e±{k) — Uq ^ \t\\/ 3 + 2[cos(/c • Si) + cos(/c • ^2) + cos(/c • ^3)] 



(4.5) 



which are nondegenerate except at the zone corners. The lowest eigenvalue occurs in s+{k) 
at k = 0, 

Eo{T) = uo- 3\t\ , (4.6) 

and dominates the partition sum (4.1) in the limit L — 00. The defect unbinding tem- 
perature Ta is determined by the condition Eo{Ta) = 0. To proceed further, we need to 
determine uq and t. 

We assume the tunnelling process between sites of the honeycomb lattice can be mod- 
eled by a one-dimensional quartic potential in the coordinate x connecting two honeycomb 
lattice sites. 



V{x) 



262 



x-lb 



x+lb 



(4.7) 



where the potential vanishes at the lattice sites. The frequency co is fixed by equating the 
maximum to Ae > 0, the energy difference between the edge and centered interstitials in 
Table 2, 



ed b 



(4.8) 
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A standard quantum mechanical calculation [43] then gives 



E± = Tuj 



2 V 27r V T 



1 



(4.9) 



for the two lowest lying levels, with 



6 



(4.10) 



The WKB exponential factor Sq is just the "kink energy" of a defect with stiffness which 
moves between the two sites as a function of z in the path integral (1.2). Upon identifying 
the splitting in (4.9) with t in Eq. (4.3), we find 



(4.11a) 



To a zeroeth approximation, we have uq = Ed, the energy of the centered interstitial 
computed in Sec. 3. In principle, uq should be corrected by the zero point energy of a two- 
dimensional quantum oscillator, i.e., by twice the first term of Eq. (4.9). This represents 
the entropy of harmonic fiuctuations of the defect. There is also, however, a negative 
contribution of this form from the entropy of the flux line lattice itself [25] which should 
have approximately the same magnitude. We assume for simplicity here that these two 
contributions simply cancel. 

It is convenient now to set 

Uq = aAe , (4.11b) 

where according to Table 3, a ^ 70. Upon using Eqs. (4.11) and (4.10), we find that the 
condition Eo{Td) = takes the form 




(4.12) 



which is solved by So/Ta = 0.085. The assumption Sd ^ Si, then leads to 

Td = 0.30y/^ao , (4.13) 

i.e., Eq. (2.24) with C3 = 0.30. The substitution cl = 0.15 - 0.30 in Eq. (2.11), however, 
shows that the melting temperature is significantly (about an order of magnitude) 
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smaller than our estimate of in the regime Bd < B < B^- Evaluation of the high 
field formula Eq. (2.26) at B = B^ is only slightly more encouraging: Using Table 3 for 
the centered interstitial energy leads to T^, = 0.079 Sodo, while the Lindemann criteria at 
B = B X gives = 0.02 — 0.09 £0*^0 ~ [44]. We conclude that the supersolid probably 
only exists for B > Bx , as indicated in Fig. 5. 

Of course, our estimate for So/Td is so small that it casts doubt on the validity of 
the tight-binding model and the WKB approximation, which are only strictly correct at 
low temperatures. More accurate band structure estimates of for B < B^ (including 
B < Bel) would certainly be of interest. 

4.2 Nature of the Transition at 

As discussed above, a transition to a supersolid below the equilibrium melting tem- 
perature of the flux crystal becomes possible when B > Bx ■ Once the defects proliferate, 
interactions between them become important, as in the closely related problem of vortex 
penetration just above Hd. To model this process, we use a continuum coherent state path 
integral representation of the partition function, similar to one used for flux lines near Hci 
[7]. The details of the lattice of preferred sites are neglected, although these could easily 
be taken into account if necessary. The grand canonical partition function describing the 
interstitial degrees of freedom then reads 

Zgr - J Vi^iVi^^e-^^^^"^*^ (4.14) 

with 

' " ^ ^ ^ ~ V (4.15) 



S[^,^*] = J (fr dz 



2£d 



Here 'i/ji{f, z) and il^lif^z) represent interstitial line creation and annihilation operators, 
and the r oc (T^ — T) is the defect chemical potential. The areal density of interstitials rii 
is given by 

ni = |V'i(r,^)|2 , (4.16) 
and the couplings u and v represent the effect of interactions. 

The nature of the transition to the supersolid which occurs with decreasing r depends 
crucially on the sign of the quartic coupling in Eq. (4.15). If tt > 0, then a continuous tran- 
sition results, and rid ^ (T — Td) up to logarithmic corrections. Both the coefficient and the 
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logarithmic correction can be calculated as in [7]. If tt < 0, the transition to the supersolid 
becomes first order. A first-order transition is possible because the microscopic two-body 
interaction between centered interstitials (represented by u\'ip{f, z)]"^ in Eq. (4.15)) is at- 
tractive, at least at large distances; see Section 3. The unbinding temperature < 
is probably too low to allow for a significant thermal renormalization of u, which could 
in principle be driven positive by entropic effects. The first-order transition described by 
(4.15) with u negative is discussed in Ref. 45. 

A finite density of proliferating defects is detectable, at least in principle, via a neutron 
diffraction experiment which precisely determines the temperature dependent magnitude 
of the six smallest reciprocal lattice vector {G{T)} of the Abrikosov fiux array for fixed 

— * — * 

magnetic field B. Note that B will also be, in general, temperature-dependent for fixed 

— * 

external field H. The number of unit cells per unit area ric is related to G according to 

nc{T) = V3G^{T)/87T^ . (4.17) 

If vacancies or interstitials only exist in small closed loop pairs, as in Fig. lb, ric must be 
exactly equal to the areal density of vortices, tt-c = uq = B/(f)Q. Above the proliferation 
temperature T^, interstitials dominate over vacancies, and we have 

= no-nc(T)>0 . (4.18) 

A plot of nc{T)/no is thus a direct measure of the density of defects, with nc{T)/no = 1 
for T < Trf, and nc{T)/no < 1 when T > T^. 

A nonzero defect order parameter {ipiir.z)) necessarily implies that the boson or- 
der parameter (1.6) is nonvanishing, because wandering vacancies or interstitial defects 
catalyze enhanced entanglement of the underlying vortex crystal. The vortices are thus 
simultaneously crystalline and entangled when < T < T^, as discussed in the Introduc- 
tion. Once the equilibrium concentration of one type of defect is nonzero, all defects will 
proliferate in at least small concentrations. Consider, in particular, the process shown in 
Fig. lb, the creation of vacancy-interstitial pairs. This can be modeled by adding terms 
to the action (4.15) as follows: 



S d^rdz 
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(4.19) 



The first two terms are the vacancy propagator and chemical potential, while the last 
allows pair creation and destruction with probability proportional to g [46]. Because the 
vacancies are unfavorable relative to interstitials, ry will remain positive just above T^. 
We see, however, that a nonzero {'i{ji{r, z)) acts like an ordering field on 'il;y{r,z), i.e., 
(V't; (r, 2) ) 7^ for T > Td- In more physical terms, vacancies can appear because their 
unfavorable energy in isolation is compensated by nearby thermally created interstitials. 

We also mention an exotic type of supersolid which is possible at least in principle: 
Suppose that the split interstitial, or some similar defect had the lowest energy. The three 
distinct orientations of this defect within its hexagonal cell represent an internal degree 
of freedom for the corresponding "boson" world lines. Should such "ribbon-like" defects 
proliferate in the solid, the resulting fluid of lines would be a "quantum rotator" liquid, 
with the same potential broken symmetries as a quantum three-state Potts model when 
interline interactions are taken into account. 

4.3 Nature of the Transition at 

It is important to distinguish between Type I and Type II melting into a liquid phase, 
according to whether vacancy and interstitial defects have already proliferated — see Fig. 4. 
Consider the standard Landau argument for a first-order transition starting from the flirx 
liquid state: Provided fluctuations suppress crystallization well below the mean fleld Hc2, 
one can expand the local BCS condensate density in Fourier components of the incipient 
density wave, 

(|^BCs(r,^)P)=Po + X;Pc5e^«- , (4.20) 

G 

where the {G} are reciprocal lattice vectors in a plane perpendicular to z. The free-energy 
difference between the liquid and crystalline phases can then be expressed as a Taylor 
series in the order parameters {Pq}, 

1 ^ 

^^=^rGY.\PGf + w Pg.Pg,Pg, + ---^ (4-21) 

Gi+Gj+Gk=0 

where we have included only the flrst ring of six smallest G's. The crucial element is the 
third-order term allowed by the symmetry of a triangular lattice, which leads to a flrst- 
order transition within this mean fleld theory [2]. The magnitude of the smallest G's in 
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Eq. (4.20) would be completely determined by the magnetic field for Type I melting, but 
would depend on the incipient density of vacancies or interstitials in the Type II CctSGj ctS 
discussed in Sec. 4.2. 

Tesanovic [42] has suggested that when the freezing transition approaches Hc2i as it 
will at low temperatures or if the fluctuations are weak, the transition will be to a charge 
density wave state with a wavevector not simply related to the magnetic field. Such a state 
is similar to the "supersolid" phase discussed here, and presumably evolves continuously 
into it as the melting field falls well below Hc2- The line of supersolid phase transitions 
may also be related to the solid- phase transition suggested by Glazman and Koshelev [24] , 
above which phase coherence is lost in the z direction. See Sec. 4.4. 

Note that generically, if a melting transition is weakly first order, it is likely to be 
to an incommensurate solid phase as the wave-vector-dependent tq in Eq. (4.21) will in 
general not have its minimum at any particularly simple wavevector. 

The distinction between Type I and Type II behavior also affects dislocation mediated 
melting theories, which start from the ordered phase just below T^. Reference 8, for 
example, studies effects of dislocation loops which are confined to the plane spanned by 
their Burger's vector and the 2;-axis. This restriction implicitly assumes Type I melting, i.e., 
that no vacancies or interstitials are present at long wavelengths. The absence of vacancy 
or interstitial lines means that only "glide" motions of dislocation lines are allowed along 
the time-like z axis, which is equivalent to a planarity restriction for the three-dimensional 
vortex loops. A small concentration of proliferating vacancy or interstitial lines would 
allow "climb-like" distortions of a vortex loop, as the loop absorbs or emits these defects. 
Type II dislocation mediated melting would thus involve arbitrary nonplanar dislocation 
loop configurations, in contrast to the planar loops associated with Type I melting. 

4.4 Transport Properties 

Finally, we compare briefly the resistive properties of the supersolid, crystal and vortex 
fluid phases. 

In the absence of pinning by random inpurities, all the above phases dissipate currents 
perpendicular to the macroscopic magnetic field since the vortex lines can move freely 
provided — in the solid-the motion is uniform. In contrast, the vortex crystal is a linear 
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superconductor for currents parallel to the z direction. Concommitantly, it can screen 
additional magnetic fields normal to the z direction. This will not occur in a semi-infinite 
system with a planar surface since the magnetic fields can just rotate. But in a cylinder 
with the vortex parallel to the axis, a small additional azimuthal magnetic field will be 
expelled, decaying exponentially in a layer with thickness given by an effective penetration 
length [48]. 

Even in the crystal phase, jinile currents parallel to z will be dissipated by nonlinear 
nucleation of vorticity normal to z\ the details of this process in bulk samples have not, to 
our knowledge, been analyzed. 

It is interesting to note that while the vortex crystal is very anisotropic in response 
to uniform currents, it appears much more isotropic in its linear confinement of magnetic 
monopoles as discussed in the Introduction. 

Both the vortex supersolid and vortex fluid phases will respond like normal metals to 
currents parallel to although there will be additional nonlinear dispersion and nonlocal 
effects [49] . The linear resistivity may also be extremely anisotropic if the vortex lines are 
rather straight. In the supersolid phase, the dissipation will be dominated by the fluid of 
interstitial (or vacancy) lines while the underlying lattice will not move in response to a 
small current in the z direction. 

Pinning by random impurities, oxygen vacancies or other defects strongly affects the 
resistive properties of the mixed state of superconductors. If the pinning is very weak, the 
thermodynamic properties will be little affected but the long-range translational order of 
both the crystal and supersolid phases will be destroyed on long distances, resulting in 
a large but finite positional correlation length. The crystal phase may be replaced in its 
entirety by a truly superconducting vortex glass phase [12] with vanishing linear resistivity. 
In contrast, the vortex fluid phase is not much affected by weak pinning. Because the 
defects are fluid-like in the supersolid phase, they also will not be strongly affected by 
weak pinning, provided they are sufficiently dense [50]. The resulting phase will then still 
be metallic with nonzero resistivity. Thus, it appears likely that the vortex-glass phase 
transition for weak pinning should occur at what was the defect proliferation transition, T^^ 
in the pure system. If the vortex lattice melting is Type /, then the putative vortex-glass 
transition should probably occur at the of the pure system. 
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The nonlinear response of the weakly pinned supersolid will involve motion of both 
the underlying lattice, including plastic flow involving dislocations, as well as vacancies 
and interstitials. Thus, even if the defect lines themselves are pinned, and the supersolid is 
in a vortex glass phase, the nonlinear response might distinguish it from a weakly pinned 
crystal. 

In the presence of sufficiently strong pinning, both equilibrium phase transitions — 
defect proliferation and melting — will be destroyed. The vortex glass transition, if it exists 
at finite temperatures, will then become second order with the resistance vanishing con- 
tinuously as the temperature is decreased. In this case, the distinctions between the low 
temperature solid phases will disappear as the extent of any crystaline order will be very 
short range. There are, however, several caveats. Firstly, it is likely that hexatic bond- 
orientational order can persist out to much longer distances than positional order [51], and 
perhaps even to infinite distances so that distinct hexatic and vortex glass transitions could 
occur. Another possibility is that more than one type of vortex glass phase could exist, 
one with residual shear elasticity, as hypothesized by Feigel'man et al. [53] and one with- 
out shear elasticity as discussed by Fisher, Fisher and Huse [12]. Which of these phases 
actually can exist even in principle is still very unclear. This complex set of possibilities 
illustrates the subtleties once random point pinning is taken into account. 

As mentioned earlier, the basic phase boundaries may well, except for strong pinning, 
be determined by the properties of the pure vortex systems analyzed in this paper. Thus, 
a more quantitative analysis of the phase diagram in the pure system is certainly merited. 
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Appendix A: Calculation of the Dynamic Matrix in the Limit A_l — > oo. 



In this Appendix we describe the details of the calculation of the dynamic matrix. 
Starting from Eq. (3.16) the central quantity to calculate is 



I{x/X) = Y,e-'''''^'^Ko{\x{l) - x\/X) - Kq{\x\/X) , 



(Al) 



where we set A_l = A in what follows. Upon using the integral representation of the Bessel 
function Kq 

^ r°° . r ^2 1 

(A2) 



1 f 

2 Jo 

and Ewald's generalized theta-function transform [39] 



4tA2 



-iq-x{l) 



exp[—\x{l) — x\'^t] = — - > e 

G 



^e-^(«"*+^)-^exp[-|^+G|V4t] , (A3) 



where Oc is the volume of the unit cell, one can rewrite the first term in Eq. (Al) as 

hix/X) = J]e-*«"'-^«Ko(|^'(0 - x\/X^) 
I 

\x{l) - xp" 



1 r°° 
^ Jo , 



-iq-x{l) 



exp 



4tA2 



2n 



X^ / dTe-^Y.e-'^^+^>^exp[-TX^\q + G\ 



(A4) 



Now we split the r-integration by introducing an arbitrary Ewald split-parameter e 

\xil)-x\^' 



^ Jo , 
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-e(A2|Q-+G|2+l) 



1 + A2|g^|2 ' 



(A5) 



where we have already performed the r-integral in the second term. 



We are interested in the limit of large London penetration depth A. Upon choosing 
e ~ A~^ we take the limit A — > oo with eA^ staying finite. Then one gets for the second 
term in Eq. (A. 5) 

(A6) 



1 + A2|g^|2 • 
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For the first term we make the substitution y = \x{l) — x\'^/{AtX'^) and get 

\x{l) — a^p"' 



1 f^dy ^ 



4yA2 



g v^-iq-Sil) 



(A7) 



In the limit A ^ oo we get exp 



x{l) — x\ 



1 since y is bounded from below by 



o = 1^(0 ~ ^P- Hence the first term can be written in terms of the exponential 



4eA2 

integral function which is defined by 



Ei{x) = -Ei{-x) = / 

J X 



In summary one gets 



lim /i(f/A) = ly"e-"^<^^Ei (■ 

A— >oo 2 \ 



y 



x{l) — x\ 



(A8) 
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-i{q+G)-x[ 
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In order to get I{\x\/X) one has to subtract 

roo 



Ko{x/X) 



ln(a;/2A) 



(AlO) 



from /i(|x|/A) in the limit A — >^ oo. Therefrom one can now calculate Ca/3 and expand in 
powers of the wave vector q. The result is 



Stt^A^ 



.2 ■ 8 



(All) 



where we have taken l/(4eA2) = tt/oc for the Ewald separation parameter. 



Appendix B: Variational Calculation of the Vacancy Formation Energy 

In this appendix we give the details of the variational calculation for the vacancy 
formation energy. Upon using the long wavelength approximation (3.18) for C°'^{q), one 

(2) 

finds three contributions to -Ey^, the vacancy relaxation energy at constant lattice spacing, 

E^'liiuie)}) = 2eoih + h + Is) (Bl) 
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with 
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= ''dqq-'f{q) + ^j^ " dqqf{q) , (B3) 



and 
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2. ^ - - (g • fc)(g^ - 2g • + fc^) ^ 



^(g-fc-/c2)(g2_g.^) 1 



{q^ -2q-k + k^)q^k^ 8' ' g^p 

(B4) 



where /(g) = /(I + eg + lig^) and we have approximated the hexagonal Brillouin zone by 
a circle of equal area with radius - 



The third contribution can be split into Is = /3a + Izh + Izc where 

T f'" W f'\^l. 2^ r {qk cos if -k'){q^-qk cos ip) 

ha = a.^J^ qdqj^ k dk J _J<p f{p) f{k) _ 3,/, eos + " ^^'^ 

Upon using 
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one arrives at 
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All these integrals diverge at small momenta. This is due to the fact that we have taken 
the penetration depth A = Aj^ to be infinite. Upon reintroducing the lower cutofi^ A~^ the 
leading contribution (which diverges in the limit A — > oo) is found to be 



/i = ^ ln(A;dA) + 0(1) , 
h = -f In(fcdA) + C(l) , 
Is =0(l,(lnA)/A2). 



Hence one finds 



(BIO) 

(Bll) 
(B12) 
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Since we are interested in the limit A ^ oo we must choose / = 1 in order to cancel the 
ln(/cdA) divergences in the final expression for the vacancy-free energy. 

Now we have to look at the subleading terms. We set / = 1 and do variational 
calculation for c and d. One gets 
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c = 



d = 
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o.ieiv^, 
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The result for the vacancy relaxation energy at constant lattice spacing becomes 
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which corresponds to a free energy of 
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TABLE 1 

Estimates for Melting and Defect Unbinding Transitions 



Regime T^(S) Td(S) 

^ B 0.5£oc^o/87rV3 const, x sqcJo/ ln(27rS/Sx) 

Ba<:B ^ Sx ciVi^(0o/5)V2 C3v^(0o/S)V2 

B ^ Bel cIv^Ai (^) e-^(^-/B)^/^ const, x V^iX± (^) e-^^^-/^)' 



/2 
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TABLE 2 



Exact and numerically calculated energy per particle of the perfect flux line lattice for 
different system sizes. N is the number of flux lines. The energies are measured in units 

of Eo = 260 . 



-E'exact -E'num 

30 0.85583194 0.85583216 

120 1.20240554 1.20240652 

270 1.40513809 1.40514051 

480 1.54897913 1.54898290 

750 1.66055090 1.66055687 

1080 1.75171168 1.75172015 

1470 1.82878705 1.82879859 
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TABLE 3 



Defect energies at constant line density for the symmetric vacancy (V6), "squeezed" va- 
cancy (V2), edge interstitial (EI), centered interstitial (CI), split edge interstitial (SEI), 
and split centered interstitial (SCI) configurations. The results are obtained from a molec- 
ular dynamics type of calculation for different system sizes with N flux lines. The energies 
are measured in units of Eq — 2eo. 



N 


Eye 


Ey2 


-E'ei 


-E'ci 


-E'SEI 


-E'sci 


30 


0.11857 


0.11394 


0.07156 


0.07218 


0.07156 


0.07166 


120 


0.12204 


0.10892 


0.07358 


0.07274 


0.07359 


0.07295 


270 


0.12381 


0.10825 


0.07392 


0.07293 


0.07392 


0.07310 


480 


0.12416 


0.10783 


0.07403 


0.07299 


0.07402 


0.07315 


750 


0.12429 


0.10781 


0.07408 


0.07304 


0.07405 


0.07319 


1080 


0.12434 


0.10772 


0.07410 


0.07309 


0.07411 


0.07324 


1470 


0.12437 


0.10780 


0.07412 


0.07307 


0.07411 


0.07325 
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Figure Captions: 



Fig. la: Dislocation loop in a flux line solid. Dashed lines represent vortices just behind 
the plane of the figure. Such loops lie in the plane spanned by their Burgers vector and 
the z-axis. The orientation of the three triangles is the same, showing that the loop has 
only a small effect on the orientational order. 

Fig. lb: Vacancy-interstitial pair in a flux line solid. Dashed lines represent vortices just 
behind the plane of the flgure. Unlike the dislocation loop in Fig. la, this loop is not 
constrained to lie in a single plane. 

Fig. 2: Vacancy line (thick dashed curve) meandering through a vortex crystal. The full 
lines show the flux lines which are in the same plane as the meandering vacancy. The 
dashed lines represent the flux lines in the neighboring plane. 

Fig. 3: Lowest energy contribution to the order parameter correlation function on the solid 
phase (1.4), this inserts a flux head and tail into the crystalline vortex array. Dashed lines 
represent a row of vortices slightly behind the plane of the page. A vacancy is created 
at "time" 2;, propagates and is destroyed at "time" z' . The energy of the "string" defect 
connecting the head to the tail increases linearly with |rj^ — fj^ | and leads to the exponential 
decay of G{f^, r^, z, z'). Physically, this represents conflnement of the magnetic monopoles 
represented by the flux head and tail. 

Fig. 4: Two distinct scenarios for vortex crystal melting with increasing temperature. In 
type I melting a flrst order transition separates a line crystal with ^ from a flux liquid 
with t/^o 7^ 0. In type II melting, both order parameters are nonzero in the intermediate 
"super solid" phase. 

Fig. 5: Schematic phase diagram of a clean high temperature superconductor. The "super- 
solid" phase is shown as the shaded region. In the presence of random pinning, a possible 
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vortex glass transition, at which the resistance vanishes, would roughly follow the melting 
boundary for low fields and the crystal-supersolid boundary in high fields. 

Fig. 6: Initial defect configurations for a symmetric vacancy, centered interstitial, edge 
interstitial, split centered and split edge interstitial used in the numerical calculation of 
the defect formation energies. 

Fig. 7: Starting from the initial symmetric vacancy configuration, shown as squares, the 
fiux line lattice relaxes first into a symmetric configuration, shown by the diamonds. This 
configuration, however, is just a saddle point and unstable against squeezing the vacancy 
along one of the three symmetry axis. One of three degenerate stable final configurations 
is represented by the circles. 

Fig. 8: Relaxed configurations for centered (squares) and edge (circle) interstitial. The 

edge interstitial is unstable with respect to a "buckling" mode perpendicular to the edge 
of the triangle. The final stable configuration is the centered interstitial configuration. 

Fig. 9a: System size dependence of the formation energies for centered (solid line), split 
centred (dot-dashed line) and edge (dashed line) interstitials. The energy in units of Eq 
is plotted versus the total number of flux lines N. The lowest energy configuration is the 
centered interstitial. 

Fig. 9b: System size dependence of the formation energies for symmetric (solid line), and 
"crushed" (solid line) vacancy. The energy in units of Eq is plotted versus the total number 
of flux lines N. The symmetric vacancy configuration (6-fold symmetry) is unstable against 
squeezing it along one of the three symmetry axis. The final stable configuration is the 
"crushed" vacancy with only a two-fold symmetry. 

Fig. 10: Relaxed configuration for an edge (open circles) and split edge (open squares) 
interstitial. The two cofigurations differ by a shift along the edge of the triangle. Because 
of their energy difference being small gliding along the edge of an triangle is a low energy 
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process. Note, however, that both configurations are unstable to buckhng, resulting in the 
centered interstitial configuration. 

Fig. 11a: Interaction energy of two centered interstitials in units of Eq versus the distance 
measured in units of the lattice constant ao for two different directions, (i) The vector 
connecting the two centered interstitials points along a, one of the basis vectors of the 
unit cell (solid line), (ii) The vector connecting the two centered interstitials points along 
b, which is perpendicular to a (dashed line). Whereas the interaction in the a-direction 
is attractive at small distances it is repulsive in the perpendicular b-direction. In the 
b-direction there is a minimum in the interaction energy at a distance of dmin ~ 2.5a. 

Fig. lib: Interaction energy of two "crushed" vacancies in units of Eq versus the distance 
measured in units of the lattice constant a. The vector connecting the vacancies is pointing 
a) along the two-fold symmetry axis of the vacancy (solid line) and b) perpendicular to it 
(dashed line). 

Fig.l2: Basis vectors {6i} for the honeycomb lattice of centered interstitial sites with 
lattice constant b. The A and B sublattices are indicated by open and closed circles. 
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